Thermal Pure Quantum States at Finite Temperature 
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An equilibrium state can be represented by a pure quantum state, which we call a thermal 
pure quantum (TPQ) state. We propose a new TPQ state and a simple method of obtaining it. 
A single realization of the TPQ state suffices for calculating all statistical-mechanical properties, 
including correlation functions and genuine thermodynamic variables, of a quantum system at finite 
temperature. 
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The possibility of extracting statistical-mechanical in- 
formation from a pure quantum state has been inten- 
sively discussed in the context of the foundation of sta- 
tistical mechanics [ll-0|- As we shall demonstrate here, 
it also has a potential significance for a new formulation 
of statistical mechanics, and for a novel calculation tech- 
nique. 

As an illustration, let us consider a closed quantum 
system composed of N spins, which is enclosed by adia- 
batic walls. In the ensemble formulation, its equilibrium 
properties are described by the microcanonical ensemble, 
which is specified by E (energy), N, and so on. The cor- 
responding subspace (energy shell) in the Hilbert space 
Hat is denoted by £e,n- Let us consider a random vec- 
tor = c^\i^) in £e,n, where {l?^)};/ is an arbitrary 
orthonormal basis set oi £e,n, denotes the sum over 
this basis, and {c^}i, is a set of random complex numbers 
drawn uniformly from the unit sphere X^i- ki'P ~ 1 t^*^ 
complex space of dimension dim £e,n- It was shown in 
Refs. that almost every such vector gives the cor- 

rect equilibrium values of a certain class of observables 
A by (iplAlip). This property was proved in Refs. 0, 0| 
for observables of a subsystem, which is much smaller 
than the whole system. The case of general observables, 
including observables of the whole system (such as the to- 
tal magnetic moment and its fluctuation), was analyzed 
in Refs. 0, It was shown that the above property 
holds not for all observables but for observables that are 
low-degreepolynomials (i.e., their degree ^ N) of local 
operators j3|. We here call such observables mechanical 
variables. We assume that all mechanical variables are 
normalized in such a way that they are dimensionless. 

For conceptual clarity, we call generally a pure quan- 
tum state that represents an equilibrium state a thermal 
pure quantum state (TPQ state). Stating more precisely 
for the case where a state \ip) has random variables (such 
as the random vector discussed above), we call {ip) a TPQ 
state if for an arbitrary positive number e 
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as — > cxD. The above inequality means that for large 
A^ getting a single realization of a TPQ state is sufh- 
cient, with high probability, for evaluating equilibrium 
values of mechanical variables. The vector c^lv) of 
Refs. [ll-Ql is a TPQ state. However, important problems 
remain to be solved. Most crucially, genuine thermody- 
namic variables, such as the entropy and temperature, 
cannot be calculated as (V'l^lV') because they are not 
mechanical variables Moreover, one needs to prepare 
a basis of £e,n to construct J^l <^i^W)- Since this 

is a hard task, such a TPQ state is hard to obtain. 

In this Letter, we resolve these problems by propos- 
ing a new TPQ state, a novel method of constructing it, 
and new formulas for obtaining genuine thermodynamic 
variables. This novel formulation of statistical mechan- 
ics enables one to calculate all variables of statistical- 
mechanical interest at finite temperature, from only a sin- 
gle realization of the TPQ state. We also show that this 
formulation is very useful for practical calculations. 

New TPQ state - We consider a discrete quantum sys- 
tem composed of A^ sites, which is described by a Hilbert 
space TLn of dimension D = A^, where A is a constant 
of 0(1). [For a spin-1/2 system, A = 2.] Our primary 
purpose is to obtain results in the thermodynamic limit: 
A^ — >■ oo while E/N is fixed. Therefore, we hereafter 
use quantities per site, h = H/N (where H denotes the 
Hamiltonian), u = E/N, and {u;N) instead of {E,N). 
[We do not write explicitly variables other than u and 
A^, such as a magnetic field.] We assume that the system 
is consistent with thermodynamics in the sense that the 
density of states g{u; N) behaves as [H] 



-l{u; N) = exp[Ars(u; A^)], /3'(w; N) < 0. 
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for every mechanical variable A. Here, P(a;) denotes the 
probability of event x, {-^e n denotes the ensemble aver- 
age, and ri^{N) is a function (of A^ and e) which vanishes 



Here, s{u; N) is the entropy density, which converges to 
the A^-independent one s{u; oo) as A^ — > oo, /3(u; A^) = 
ds{u; N)/du is the inverse temperature, and /?' = 8(3/ du. 
These conditions are satisfied, for example, by spin mod- 
els and the Hubbard model. Since D is finite, j3 may be 
positive and negative in lower- and higher-energy regions, 
respectively. We here consider the former region. 

We propose the following TPQ state and the proce- 
dure for constructing it. First, take a random vector 
■f/^o) = X^i^il*) from the whole Hilbert space TLn- Here, 
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is an arbitrary orthonormal basis of 'HN^ and {ci\i 
is a set of random complex numbers drawn uniformly 
from the unit sphere ~ 1 of the £)-dimensional 

complex space. Note that this construction of random 
vectors is independent of the choice of the orthonormal 
basis One can therefore use a trivial basis such 

as a set of product states. Hence, |V'o) can be generated 
easily. On ther whole, the amplitude is almost equally 
distributed over all the energy eigenstates in this state 
(as is easily seen by choosing the eigenstates of h as the 
basis Thus, the distribution of energy in |'0o) is 

proportional to g{u\ N). We wish to modify this distribu- 
tion into another distribution rk{u] N) which has a peak 
at an desired energy. This is easily done by operating a 
suitable polynomial of h onto jV'o} as we shall see below. 
[Operating h onto a vector is much easier than diago- 
nalizing h.] We denote the minimum and the maximum 
eigenvalues of h by Cmin and Cmax, respectively. Take a 
constant / of 0(1) such that I > emax- Starting from 
\ipo), calculate 



li'k+i) 



{ipk\h\-)pk), 

{l-h)\^Jk)/\\{l 



h)\il^k)\ 
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iteratively for k — 0, 1,2,---. From Eq. ([7]) below, uq 
corresponds to /3 = 0, i.e., g{u; N) takes the maximum 
at M = mq. We will also show that Uk decreases gradually 
down to Cniin as k is increased, i.e., uq > ui > ■ ■ ■ > 
Cmin- One may terminate the iteration when Uk gets low 
enough for one's purpose. We denote k at this point by 
^tcrm- We will show that fcterm = 0{N) at finite tem- 
perature, and that the states \ipo), ■ ■ • , IV'fetorm) be- 
come a series of TPQ states corresponding to various en- 
ergy densities, Mo,iti, ■ • ■ ,Ufct„„„. Hence, the equilibrium 
value of an arbitrary mechanical variable A is obtained 
as {ipk\A\ipk), as a function of Uk- For each realization of 
{ci}i, a series of realizations of TPQ states is obtained. 
We will show that the dependence of {tpk\A\'^k) on {ci}i 
is exponentially small in size N as N increases. There- 
fore, only a single realization suffices for getting a fairly 
accurate value. When better accuracy is required, one 
can take the average over many realizations. 

We now show that the states obtained with the above 
procedure are TPQ states. Since |V'o) is independent 
of the choice of the basis, we take the set of energy 
eigenstates {|n)}„ as {\i)}i in order to see properties 
of \'4>k) (although we never use such a basis in practi- 
cal calculations). After /c-times multiplication oi I — h, 



IVfc) cx il-h)''\^o) 



]cn{l ~ en)''\n), 
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where h\n) — e„|n). Let us examine how the energy 
density u distributes in this state. The (unnormal- 
ized) distribution function of u is given by rk{u; N) = 



(5~-^Xjj"|c„p(Z — en)^*^, where Sr = o(l) and the sum 
is taken over n such that e„ lies in a small interval 
[u — Sr/2, u + (5r/2). Since the density of states g{u; N) 
is exponentially large in size N, rk{u; N) converges (in 
probability) exponentially fast to its average. Hence, 



rk{u;N)^D-^exp[NUu;N)], 



(6) 



where ^^.{u; N) = s{u; N) + 2Kln{l - u) with k = k/N. 
Hereafter we often denote k dependence by k, e.g., we 
express Uk as u^. Note that ^k{u;N) does not depend 
on {ci}i, because the dependence vanishes when we have 
dropped negligible terms in Eq. £^^{u; N) takes the 
maximum at u* which satisfies 



Piu:;N) = 2K/il-u:). 



(7) 



Since /3(u*;iV) and I - m* are 0(1), we find k = 0(1), 
and hence k — 0{N). Expanding £^k.{u;N) around u*, 
and noticing 

C ^ d'Udu' = 13' {ul; N) - 2n/{l u*^ < 

from Eq. Q, we get i^{u:N) = i.{ul-N) - |e;,'|(u - 
<?l'2 + - +■■■■ Here, ^ = d^Udu^ = 

/3"(m*; iV)-4K/(Z-u*)^. Hence, rk{u; N) behaves almost 
as the Gaussian distribution, peaking at m = m* , with the 
vanishingly small variance l/iV|^"|. Let us introduce the 
density operator pk = {l — Kf'^/Tr{l~ h)'^^, which has the 
same energy distribution rk{u;N). In the ensemble for- 
mulation, pfc represents the equilibrium state specified by 
(m„;7V) because ri^{u;N) has a sharp peak. We call the 
ensemble corresponding to pk the smooth microcanonical 
ensemble (because the energy distribution is smooth). In 
a way similar to those of Refs. [1, we can show that 
for an arbitrary positive number e 



>.)<^fei^^, (8) 



2rfc«;iV) 



(^fc|^|i/;fc) =Tr[pfcA] 



(9) 



for every mechanical variable A. Here, || • || de- 
notes the operator norm 0], and the overline repre- 
sents the random average. With increasing N, 
grows at most as a low-degree polynomial of N , whereas 
'i'k(emin',N)/rk{u*^]N) decreases exponentially at finite 
temperature (i.e., for it* > emin)- Therefore, li/j/t) is a 
TPQ state for the smooth microcanonical ensemble. 

Genuine thermodynamic variables - One might think 
it impossible to obtain genuine thermodynamic variables 
like the temperature and entropy by only manipulat- 
ing pure quantum states. However, our new TPQ state 
makes it possible. In fact, by substituting for u* in 
Eq. ([7]), and using Eq. (|lip below, we obtain 



f3{u^-N)^2K/(l-u,) + 0{l/N). 



(10) 



This gives /3(uk;; N)^ with an error of 0(1/A^), as a func- 
tion of [because k and I are known parameters] . That 
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is, one obtains the temperature of the equihbrium state 
specified by (mk;7V) just by calculating u„ with Eq. ([3]). 

We can also obtain formulas with less errors. For ex- 
ample, using Eq. ^ and the expansion of ^„(u;7V), we 
have 

< = < + 0{1/N'), ul^u,- C/2mf. (11) 
Substituting u* for u* in Eq. ([7]), we get a better formula 
/?«;iV) = 2«/G-<) + 0(l/iV2). (12) 

One can evaluate ^" and easily by calculating {tpk \ {h— 
u.y\A) = 1/A^ICI + 0{1/N^) and Oj^lih - u^f\A) = 
C/N^\C^\^ + 0{1/N^). Hence, using formula one 
obtains P{u;N) (for u = Uq,u',---) with an error of 
0(l/iV2). In a similar manner, we can obtain formulas 
whose errors are of even higher order of 1/iV. 

However, /3(u; N) is the inverse temperature of a finite 
system, whereas we are most interested in its thermody- 
namic limit /3(u; oo). In general, the difference |/3(u; N) — 
/3(u; oo)| decays not so quickly as 0{1/N^). To obtain an 
even better formula for f3(u; oo), we consider C identical 
copies of the iV-site system. We denote quantities of this 
CTV-site system by tilde, such as \ipo) = \ipo)'^'^ . The 
state IV'fc) is given by oc (l — h)'^'''\'ipo), where h = 
i»(c-i) + i^ij ^ io(c-2) _^ . . . _^ ^ij) /CN. 

In the limit of C — > oo, ua approaches the canonical 
average of u in a single copy with inverse temperature 
l3{ufi; oo). At the point where l3{ufi;oo) = /3(it*;iV) is 
satisfied, we can estimate this canonical average, which 
is denoted by u^, in the same manner as Eq. (11). Then, 
we get u'^ — + 0(l/iV^), where 



2N[^- + 2k/{1 - 



We thus find 



/3(u:;oo) = 2^/(/-<) + 0(l/iV2), 



(13) 



(14) 



which gives the inverse temperature /3(m; oo) (for u = 
Mq, u", • • •) of an infinite system composed of an infinite 
number of iV-site systems. We expect that /3(u; oo) is 
much closer to /3(u; oo) than /3(u;iV), because informa- 
tion of f (u; N) in the whole spectrum range of u is in- 
cluded in /3(u; oo). [By contrast, only the information at 
the peak of ^(u; N) is included in (3{u; N).] This will be 
confirmed later by numerical computation. 

We can also obtain the entropy density s as a function 
of u and h^, by integrating (3 over u and /3mz over hz. 
For example, for an arbitrarily fixed value of hz , we have 

9-1 ^ 
S{u2p)-S{u2q) = Xl"("2£>^2^+l'"2£+2)+0(-r72)- (15) 



by generalizing Simpson's rule. Here, u stands for {u; N) 
or (7i;oo), p and q are integers, and v(x,y,z) = (x — 
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FIG. 1. Magnetization plotted against a magnetic field for 
J = —1. Solid lines represent exact results for N ^ oo, 
for various values of the energy density u Q. Circles denote 
results obtained with our formulation for N = 24. Results for 
N = 4-20 are also shown for u — — 0.3J. 



z){/3(x) -I- /3(z)}/2 -{x- z)2[a;{/3(z) - /3(y)} -|- y{/3(x) - 
I3{z)} + z{/3{y) — /3{x)}]/6{x ~y){y ~ z). We have also de- 
veloped another method of obtaining s, in which g{u; N) 
is directly evaluated from the inner products among dif- 
ferent realizations of a TPQ state 8]. 

To sum up, one can obtain a series of TPQ states and 
values of all variables of statistical-mechanical interest, 
by preparing a random vector and simply applying (/ — h) 
iteratively. That is, we have established a new formula- 
tion of statistical mechanics, whose fundamental formu- 
las are Eqs. ^ and (ITU)) . 

Numerical results - Our formulation is easily imple- 
mented as a method of numerical computation. We ap- 
ply it to the one-dimensional Heisenberg model in order 
to confirm the validity of the formulation. We take H — 
i Yld=l[^i.'^)^^i.'^+'^)-hz<yz{^)], where J = -1 (ferromag- 
netic) or -1-1 (antiferromagnetic). For N ^ oo, the exact 
results at finite temperature (i.e., u > emi„) have been 



derived for magnetization mz 



TV- 



all values of u and h 
-1 v-Af 



N 



at 



and for the correlation function 



(f){j) EE N ^ y^i(crz(j)a2(i -h j))°^^ at all values of u 
with hz — [lO|. They are plotted in Figs. [T] and by 
solid lines, where different colors correspond to different 
values of u. We calculate the corresponding results using 
our formulation, by performing numerical computation. 
The results for N = 24 are plotted by circles, where each 
circle is obtained from a single realization of TPQ state. 
According to Eq. (jS]), choice of the initial random num- 
bers {ci}i has only an exponentially small effect on the 
results at finite temperature. We have confirmed this 
fact by observing that the standard deviation, computed 
from ten realizations of a TPQ state for each data point, 
is smaller than the radius of the circles of these figures. 

Results for other values of N are plotted in Fig. [1] for 
u — —0.3 J, and in the left insets of Fig.[2]for u = —0.36 J 
at j = 2. It is seen that the TV-dependence becomes 
fairly weak for A'' > 20, and that the results for A^ = 24 
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FIG. 2. Correlation function plotted against j for J = 
+1 and hz = 0. Solid lines represent exact results for A'^ — )■ 
oo, for various values of u [lol ]. Circles denote results of our 
formulation for = 24. (Left Inset) Results for = 16-24 at 
j — 2 for u — —0.36 J. (Right Inset) at finite h^, obtained 
from a single realization of the TPQ state at T ~ 0.45J. 




u 



FIG. 3. Temperature T plotted against u for J = +1. 
Solid lines re pre sent exact results for N ^ oo, for various 
values of UM- Triangles and squares denote l/fi{u\N) 
(triangles) and l//3(u;oo) (squares) for A'' = 24, obtained 
with our formulation. (Inset) l//3(u; oo) for A'^ = 8-24. 



agree well with the exact results. As illustrated by this 
example, N should be increased in our method until the 
variation of the results with increasing N becomes less 
than the required accuracy. 

We have also computed (f){i) at finite and T, for 
which exact results are unknown. The results a.t T ^ 

0. 45 J are plotted in the right inset of Fig. [5] 

For genuine thermodynamic variables, the exact result 
for 1//3(m; oo) [lH is plotted by solid lines in Fig.[3l Cor- 
responding results for l//3{u; N) and 1//3(m; oo), obtained 
with our method with N = 24, are plotted by triangles 
and squares respectively, where each point is obtained 
from a single reahzation of the TPQ state. [We have 
confirmed again that dependence on the choice of {ci}i 
is negligibly small.] Not only f3{u; N) but also f3{u; oo) 
depend on N. However, the dependence of /3(w; oo) be- 
comes fairly weak for N > 20, as shown in the inset. 
P{u;oo) for N = 24 agrees well with the exact result, 
whereas (3{u; N) differs significantly from them for this 
value of N. We have thus confirmed that (3{u; oo) is much 
closer to oo) than /3{u;N), for finite N. Note how- 
ever that N) gives almost correct result for /3 of a 
finite system, as seen from Eq. (I12p . 

We have obtained a series of TPQ states at the discrete 
points Mo, wi, U2, • • • , Utcrm- The discrete points are dense 
enough because their intervals are 0{1/N), vanishing as 
iV — > oo. The intervals also depend on the parameter 

1. When smaller I is taken, k gets smaller to reach the 
same u and temperature T, as seen from Eq. ([7]). Hence, 
to obtain results at low T, I ~ emax is appropriate to 
reduce the amount of computation. At high T, however, 
Uk moves quickly as k increases, for such a small I. Hence, 
to obtain results at many values of u at high T, I should 
be taken larger. When computing the data for Figs. [T] 
and m we have taken I ~ Cmax- Since the values of u 
which are specified in these figures are not necessarily 



found among Ufc's, we have slightly tuned / in such a 
way that Uk can be found within O.OOIJ of the specified 
values. In these figures, and 0(j) at such u^s are 
plotted. When computing the data for Fig. [31 we have 
performed computations with two values of Z; I — Cmax 
and 5J. Both results agree well with each other. For 
better visualization, we have plotted the results with I — 
5 J (orange and purple) and those with I = Cmax (red and 
blue) in the high- and low-T regions, respectively. 

Advantages - We now discuss advantages of our formu- 
lation when used as a method of numerical calculation. 
At finite T, an exponentially large number of states are 
included in £e,n- This makes computation of eigenstates 
in Se.n pretty hard. In contrast, our method takes full 
advantage of such a huge number of states, as seen, e.g., 
in the derivation of Eq. ([6]) . As a result, using just a single 
realization of TPQ state, one can calculate all quantities 
of statistical-mechanical interest at finite T, on the solid 
theoretical basis that is developed in this Letter. More- 
over, our method is applicable to systems of any spatial 
dimensions, and to frustrated or fermion systems as well. 
Furthermore, our method costs much less computational 
resources than the numerical diagonalization. For ex- 
ample, the number of non-vanishing elements of H of 
the Heisenberg model is 0(iV2^). Since k = 0{N), the 
computational time is 0{N'^2^) in our method, which 
is exponentially shorter than that of diagonalization. In 
fact, it took only two hours to compute all data in Fig. [3] 
on a PC. Computations can be made even faster by par- 
allelizing the algorithm, which is quite easy and efficient 
because our method consists only of matrix multiplica- 
tions. 

Furthermore, our method is effective over a wide range 
of T because the rhs of Eq. (jS]) is exponentially small as 
long as s (and hence T) is finite of 0(1). In fact, Figs.[l][3] 
show that our results agree well with the rigorous results 
in a wide range of T, from T J toT ^ J. In practical 
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computations with finite N, T (= l//3(u*; cxo)) can be 
lowered as long as r kiemin', N) /rk{u*^] N) ^ 1. 

We note that the quantum Monte Carlo method may 
be much faster. However, it suffers from the sign problem 
in frustrated systems and fermion systems. The density- 
matrix renormalization group method has been extended 
to finite temperature, and the state obtained in Ref. [15 1 



might be close to TPQ states. However, its effectiveness 
in two- or more-dimensional systems is not clear yet. The 
states obtained with the microcanonical Lanczos method 
[l6| . which tried to obtain not TPQ states but eigen- 
states, might also be close to TPQ states. However, the 
method costs more computational time in Ref. [l6| than 
ours, and a method of computing T or s seems more dif- 
ficult than ours. We therefore expect that our method 
will make it possible to analyze systems which could not 
be analyzed with other methods. 

Concluding remarks - We conclude this Letter by 
making several remarks. First, one can evaluate the mag- 
netic susceptibility {drriz/ dhz)u from Fig. [l]or[2l One can 
also obtain {dmz/dhz)T with the help of Fig. [3] 

Second, \i/jk) remains to be a TPQ state after 
time evolution, since Eq. ([5]) shows that e^*/*''|-!/)fc) oc 
■^^^-ie„t/h^^f^l _ g^)fc|j2)j which is just another realiza- 
tion of lipk)- 

Third, our formulation is advantageous to analyses 
of phase transitions. As an example, consider the case 
where the energy density u(T; N) for iV — ^ oo is discon- 
tinuous at the transition temperature Ttr of a first-order 
transition. Then, the specific heat c = du/dT diverges at 
T = Ttr- If one used the canonical formalism, where T is 
an independent variable, calculation of u{T; oo) would be 
hard around T = Ttr- In our formulation, by contrast, u 
is taken as an independent variable, and c is obtained as 
c= — /?^/(9/3/9u)jv from /?(m;oo). The function /3{u;oo) 
is continuous even at the transition point, where it takes 
a constant value 1 /Ttr in a finite interval of u correspond- 
ing to the phase coexistent region Hence, f3{u;oo) 
can be calculated more easily than u(T; oo). In fact, one 
can identify a first-order transition by simply observing 
that the rhs of Eq. ([T4|) takes a constant value, apart from 
small deviation of 0{1/N'^), for multiple values of m* and 
K. Regarding a continuous transition, it can be identified 
from a singularity in c, or an order parameter to, and so 
on. One can calculate m by adding a symmetry-breaking 
field / to the Hamiltonian, and thereby computing m at 
/ — ±1/1 for small |/|. Or alternatively, without intro- 
ducing /, one can perform the 'pure-state decomposition' 
(i.e., decomposition into macroscopically definite states) 
by a pply ing the variance- covariance matrix method of 
Ref. |l3 to a TPQ state. 

Fourth, Eq. ^ can be generalized as \4>) oc Q{h)\'ipo), 
which defines other new TPQ states. Here, Q{u) is any 
differentiable real function such that Q{u)^g{u; N) has 
a sharp peak, whose width vanishes as iV — > oo, and 



Q{u)'^g{u; N) outside the peak decays quickly. Using this 
\ip), one can calculate various quantities as we have done 
using lipk)- For instance, the formula corresponding to 
Eq. is given by l3{u*-N) + 2Q'{u*)/NQ{u*) = 0. 

Finally, although a TPQ state (such as \tpk)) and the 
mixed state (such as pk) of the corresponding ensemble 
are identical with respect to mechanical variables, they 
are completely different with respect to entanglement. At 
T ^ J, for example, pk has only small entanglement 
(because it is close to the completely mixed state (1/Z?)1, 
which has no entanglement), whereas we can show that 
\ipk) has exponentially large entanglement (as previously 
shown for T oo in Ref. [l3l)- It is thus seen that an 
equilibrium state can be represented either by a TPQ 
state with huge entanglement or by a mixed state with 
much less entanglement. Their difference can be detected 
only by high-order polynomials of local operators, which 
are not of statistical-mechanical interest |3l. Il4|. 
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